Method and system for estimating directions-of-arrival in low power or low sample size scenarios

ABSTRACT

The present invention relates to a method and a system for determining the direction of arrival of one or multiple radio or acoustic waves and, more particularly, to such a method and system especially advantageous in situations where the number of available observations is small and/or the available observations are received with low power. The method significantly improves the performance of traditional subspace signal processing algorithms in the low sample size regime. The algorithm is specifically designed to provide consistent estimates even when the observation dimension has the same order of magnitude as the number of observations. This guarantees a good behavior in finite sample size situations, where the number of sensors or antennae and the number of samples have the same order of magnitude or in scenarios where the received signal power is not sufficiently high to guarantee detection via conventional one-dimensional search methods.

FIELD OF THE INVENTION

The present invention relates to a method and a system for determiningthe direction of arrival of one or multiple radio or acoustic waveformsand, more particularly, to such a method and system especiallyadvantageous in situations where the number of available observations issmall, or the available observations are received with low power.

STATE OF THE ART

In multiple applications, there is an especial interest in determiningthe direction of arrival (DoA) or also, angle of arrival, of one ormultiple radio or acoustic waveforms as seen from a particular receiver.The estimation of the impinging angles becomes feasible when an array ofmultiple sensor or antenna elements spatially distributed in aparticular region is employed. Each of the sensors is equipped with aradio-frequency/acoustic front-end and its corresponding signalprocessing technology. The waveforms or signals received by each of thearray elements are processed together by a central processing unit thatdetermines the direction of arrival of the different sources in thescenario. The sources might be transmitters themselves, or they might bereflectors of the acoustic or radio waveforms generated by othersources.

For example, in radar applications it is useful to estimate the exactposition of different targets as seen from a given reference point. Thiscan be done by sending a particular radio-frequency waveform andestimating the delay and the direction of arrival from which the signalis received after bouncing off from different objects in the scenario.The use of an antenna array allows for the estimation of the exactposition of the target without resorting to triangularizationtechniques.

On the other hand, in some other applications such as mobilecommunications, there is an interest in estimating the position of asignal source, such a mobile station. Again, this can be effectivelydone by using an antenna array as a receiver. The direction of arrivalof the different mobile stations that are simultaneously transmittingtowards a base station can be estimated using the techniques describedherein.

There is a vast amount of methods for determining thedirection-of-arrival (DoAs) of the distinct sources, from which the mostimportant contributions are: the Capon estimator [J. Capon, “Highresolution frequency-wavenumber spectrum analysis”, Proceedings of theIEEE, vol. 57, pp. 1408-1418, August 1969], the Maximum Likelihood (ML)method [P. Stoica and A. Nehorai, “Performance study of conditional andunconditional direction-of-arrival estimation”, IEEE Transactions onASSP, vol. 38, pp. 1783-1795, October 1990], and the super-resolutiontechnique known as MUSIC (“MUltiple SIgnal Classification”) [R. Schmidt,“Multiple emitter localization and signal parameter estimation”,Proceedings of the RADC, Spectral Estimation Workshop, Rome (NY), pp.243-258, 1979. Reprinted in “Modern Spectral Analysis II”, S. B. Kesler(ed.), IEEE Press, New York, 1986].

If the number of samples is sufficiently high, the ML method gives thebest performance, although this is at the expense of a much highercomputational complexity (the ML solution is computed by optimizing acost function in multiple variables and the solution is obtained from amulti-dimensional search). Both the Capon and the MUSIC spectralestimators are much more efficient in terms of computational burden,because the solution is obtained from a one-dimensional search. Thedifference between these two methods is the fact that the MUSIC method,contrary to the Capon solution, exploits the eigenstructure of thespatial correlation matrix of the observations. In practice, it has beenshown both theoretically and experimentally, that DoA estimationtechniques based on the eigenstructure methods are superior to othertraditional one-dimensional search based methods such as the Caponestimator.

On the other hand, V. Girko introduced the theory of G-estimation [V.Girko, “An Introduction to Statistical Analysis of Random Arrays”, TheNetherlands: VSP, 1998]. This theory provides a systematic approach forderiving estimators with excellent properties in the small sample sizeregime.

The main drawback of subspace-based methods, such as MUSIC, is the factthat their capability of detecting closely spaced targets dropsdramatically when either the number of available samples or the power ofthe received waveform fall below a certain threshold. Several methodshave been proposed to alleviate the bad performance of the MUSIC methodin the low power and/or sample size region, although most of thesemethods can only be applied to some particular array geometries. Forexample, the polynomial-root MUSIC method or the state-space approachcan only operate on uniform linear array architectures (that is to say,on arrays where all the sensors/antennas are equispaced and locatedalong a straight line). On the other hand, it has been pointed outseveral times in the literature that a weighted version of the MUSICmethod might help improve the performance beyond the performancebreakdown threshold. However, up to date, there is no clear indicationon how to choose these weights in order to optimize the performance ofthe method in the low power and/or low sample size regime without losingthe good properties of the MUSIC algorithm in more generic scenarios.

SUMMARY OF THE INVENTION

It is a primary aim of the present invention to provide a method andsystem for overcoming the MUSIC performance breakdown effect, allowingfor proper detection of directions-of-arrival, hereinafter referred toas “DoAs”, of closely located targets in situations where the number ofavailable samples is low, or when the power of the received waveforms isnot sufficiently high to guarantee the detection of their DoA byone-dimension exploration methods. As we use the teachings of the theoryof G-estimation, our method will hereinafter be referred to as theGMUSIC method. Our method can be further complemented with someenhancements of the traditional MUSIC technique, such as polynomialrooting, beamspace projection, sequential estimation or spatialsmoothing.

The object of the present invention is a new method and system for DoAestimation for multi-sensor or multi-antenna array receivers. The mainparticularity of the proposed invention is the fact that it overcomesthe threshold effect of traditional eigenstructure methods, allowing forcorrect multiple DoA detection in situations where either the number ofavailable samples or the power of the received waveforms are notsufficiently high to guarantee DoA detection by conventional means basedon one-dimensional exploration of the angle of arrival.

The proposed method of DoA identification comprises the selection of thelowest or highest peaks of a cost function that depends on the angle ofexploration. This function is constructed from the eigenvalues andeigenvectors of the estimated covariance matrix of the observation.

Another advantage of our method, when compared with the well-knownimprovements of the MUSIC algorithm, is that it is valid for any type ofarray of receiving elements, such as linear array, non-linear array,uniform array, circular array, triangular array, and so on.

The proposed method for determining the direction of arrival of at leastone radio waveform, comprises the steps of: receiving at least onewaveform at M receiving elements; obtaining one detected signal at theoutput of each of the M receiving elements; processing said M detectedsignals for obtaining M processed signals; sampling at N time instantssaid M processed signals; from the samples obtained in the previousstep, creating N column vectors y(1), . . . , y(N), wherein y(n)represents, for a sampling instant n, a column vector with M rows, eachof the M rows containing the sample value associated to thecorresponding receiving element at that sampling instant n; from thecollection of N column vectors, calculating a sample spatial correlationmatrix {circumflex over (R)}, obtaining all the eigenvectors ê_(i) andall the eigenvalues {circumflex over (λ)}_(i) from the sample spatialcorrelation matrix {circumflex over (R)}; from all the eigenvectorsê_(i) and all the eigenvalues {circumflex over (λ)}_(i), obtaining a setof signal dependent parameters μ_(k); from the eigenvalues {circumflexover (λ)}_(i) and the signal dependent parameters μ_(k), constructing aset of weights; and finding the K solutions from a certain cost functionsuch that they fulfil a certain requirement.

Finding the K solutions from said certain cost function is done eitherby selecting K deepest local minima of a certain cost function, or byselecting K highest local maxima of a certain cost function.

The present invention also provides a system which comprises meansadapted for carrying out the steps of the method, as well as a computerprogram comprising computer program code means adapted to perform thesteps of the method when said program is run on a computer, a digitalsignal processor, a field-programmable gate array, anapplication-specific integrated circuit, a micro-processor, amicro-controller, or any other form of programmable hardware.

The advantages of the proposed invention will become apparent in thedescription that follows.

BRIEF DESCRIPTION OF THE DRAWINGS

To complete the description and in order to provide for a betterunderstanding of the invention, a set of drawings is provided. Saiddrawings form an integral part of the description and illustrate apreferred embodiment of the invention, which should not be interpretedas restricting the scope of the invention, but just as an example of howthe invention can be embodied. The drawings comprise the followingfigures:

FIG. 1 is a schematic representation of the receiving system utilizedfor estimating the direction of arrival of one or several waveforms inaccordance with the preferred embodiment of the present invention.

FIG. 2 is a function block diagram of the steps carried out by apreferred method of the present invention.

FIG. 3 is a graph of a particular instance of the cost function used fordirections-of-arrival “DoA” extraction as obtained in practicing thedescribed method.

FIG. 4 and FIG. 5 are graphs of the results obtained from practicing themethod or our preferred embodiment in terms of mean squared error andprobability of detection, respectively.

DESCRIPTION OF PREFERRED EMBODIMENTS OF THE INVENTION

The implementation of the present invention can be carried out accordingto the setting illustrated in FIG. 1. One or several waveforms (w_(a),w_(b), w_(c) . . . ) are received by at least one receiving element (1,2, 3, 4 . . . ). These waveforms are preferably radio-frequency ormicrowave signals. The receiving elements are preferably a sensor arraycomposed of multiple sensors or antennae. The one or several waveforms(w_(a), w_(b), w_(c) . . . ) come from respective one or severalsources, which might be transmitters themselves, or they might bereflectors of the acoustic or radio waveforms generated by othersources. The waveform or signal (d₁, d₂, d₃, d₄ . . . ) captured ordetected by each one of these receiving elements (1, 2, 3, 4 . . . ) isreceived by a pre-processing unit or front-end (30, 31, 32, 33 . . . )that pre-processes the signal (f₁, f₂, f₃, f₄ . . . ) so that it can befurther processed by the signal processing unit (40). The signalprocessing unit (40) takes samples of the signals (f₁, f₂, f₃, f₄ . . .) coming from the front-ends (30, 31, 32, 33 . . . ) at certainsynchronized time instants, and generates a collection of discrete timemultidimensional observations or snapshots [(f¹ ₁, f¹ ₂, f¹ ₃, f¹ ₄ . .. ), (f² ₁, f² ₂, f² ₃, f² ₄ . . . ), . . . , (f^(N) ₁, f^(N) ₂, f^(N)₃, f^(N) ₄ . . . )], wherein N is a natural number which represents thenumber of time instants at which samples of (f₁, f₂, f₃, f₄ . . . ) aretaken.

Assuming that M receiving elements, sensors or antennae (1, 2, . . . ,M) are available, that is to say, that the sensor or antenna array isformed by M receiving elements, at each sampling time n the signalprocessing unit (40) acquires a collection of M values (f^(n) ₁, f^(n)₂, . . . , f^(n) _(M)) coming from the M different front-ends (30, 31,32, 33 . . . ). These values (f^(n) ₁, f^(n) ₂, . . . , f^(n) _(M)) maybe real-valued or complex-valued depending on the characteristics of thereceiver. The M values are stacked into an M×1 column vector y(n), wheren is a natural number that indexes the sampling instant associated withthe observation. Hence, if the received signals or waveforms aresynchronously sampled at N different time instants, the signalprocessing unit (40) works with a collection of N column vectors orsnapshots, {y(1), . . . , y(N)}. The number N of available snapshots canbe lower than, equal to or higher than the number of sensors M.

Assuming that there are K sources, which can be transmitters, targets orreflectors of the acoustic o radio waveforms generated by other sources,in the scenario, each of these column vectors can be described asy(n)=As(n)+n(n)  (1)where:

-   -   s(n) is a K×1 column vector containing the samples of the        signals transmitted or reflected by K different sources or        reflectors.    -   n(n) is an M×1 column vector that contains the (thermal) noise        samples generated by the front-ends (30, 31, 32, 33 . . . ) at        the nth time instant.    -   A is an M×K matrix with columns a(θ₁), . . . , a (θ_(K)), where        a(θ) is the array response associated with a signal coming from        an angle θ, and θ₁, . . . , θ_(K) are the DoAs of the K        directional sources that are present in the scenario. The array        response a(θ) is an M×1 column vector, the entries of which are        functions of the angle of arrival. Since each waveform arrives        at the sensor or antenna array of M receiving elements with a        different angle of arrival θ (or direction of arrival), in the        signal represented by equation (1), each one of the K sources        contributes with a column vector M×1 a(θ₁) . . . a(θK), which is        obtained by evaluating the function a(θ) at the angles θ₁, . . .        , θ_(K).

These functions are assumed to be either known or obtained by some othermeans at the signal processing unit (40). To illustrate how this columnvector a(θ) is constructed, consider a uniform linear array of Melements and assume that a plane waveform impinges on the array from anangle of θ with respect to the broadside of the array. Assume that thesignal received at the first antenna is s(n) (that is, the transmittedsignal, disregarding any delay to simplify the exposition). Then, if thebandwidth of this signal is sufficiently low, the signal received by thesecond antenna will be s(n)exp(j2πd/λ·sin(θ)), where: j is the imaginaryunit, d is the array interelement separation, and λ is the wavelength ofthe impinging waveform. Furthermore, the signal received in the mthantenna will be s(n)exp(j2π(m−1)d/λ·sin(θ)). Hence, gathering the signalreceived at each antenna in an M×1 column vector, it turns out that onecan express this column vector as s(n)·a(θ), where again s(n) is thetransmitted signal and a(θ) is an M×1 with the mth entry equal toexp(j2π(m−1)d/λ·sin(θ)). In general, the column vector a(θ) can beconstructed for any array geometry, and it gives the array response as acomplex function of the direction of arrival θ.

As we have already said, these functions a(θ₁) . . . a(θ_(K)) areassumed to be either known or obtained by some other means at the signalprocessing unit (40).

In principle, if the receiving chains were ideal, knowing the geometryof the array of receiving elements would be enough to know a(θ).Nevertheless, in practice, the receiving chains are not ideal, but theyadd a certain amount of distortion. For this reason, in this type ofarchitectures of arrays of receiving elements, it is usual to eitherperform the step of finding out the form of the function a(θ), orcorrect the response of the receiving chains in order to obtain thedesired function a(θ). This process is sometimes called “arraycalibration”. In order to find out the form of the function a(θ),reference signals are usually introduced in the antennae or sensors,these reference signals simulating waveforms coming from the directionθ. Then, the response of the receiving chains is measured at the signalprocessing unit (40), and the result of this measures is called a(θ). Inprinciple, this process of calibration must be done for any value of θwithin the exploration margin, such that the collection of all thepossible functions a(θ) is obtained for every value of θ. However, sincethe process of calibration is done prior to the exploration itself, inpractice it is impossible to know a(θ) for any value of θ (theexperiment would have to be repeated for infinite values of θ). That iswhy, in practice, a collection of as many discrete values of θ aspossible and as close to each other as possible, is selected, in orderto evaluate and know the function a(θ) during the process ofcalibration.

Later on, while the algorithm of detection of directions-of-arrival isbeing executed, the cost function η(θ) is normally evaluated for thevalues of θ for which the function a(θ) is known, by interpolating therest of the values of θ.

As far as planar arrays are concerned, that is to say, arrays in whichthe receiving elements are distributed on a plane, the function a(θ)depends on two angles (elevation and azimuth) instead of on one angle(this being the case for linear arrays, that is to say, arrays in whichthe receiving elements are distributed on a straight line).

In summary, θ can be either a parameter or a vector or parameters, whosevalues are angles that define the direction of arrival.

The array of M receiving elements (1, 2, 3, 4 . . . ) can also becalibrated as follows: exposing a collection of reference signals to thearray of M receiving elements (1, 2, 3, 4 . . . ); measuring, at thesignal processing unit (40), the response of the received referencesignals obtained at the output of the M pre-processing units (30, 31,32, 33, . . . ), and modifying the response of the M pre-processingunits (30, 31, 32, 33, . . . ) in such a way that the signals measuredat the output of the M pre-processing units (30, 31, 32, 33, . . . ) areas close as possible to the reference signals inserted in the firststep.

FIG. 2 illustrates a preferred method for obtaining the directions ofarrival θ₁, . . . , θ_(K) from the observations {y(1), . . . , y(N)}.The proposed method uses the teachings of the theory of G-estimatorsintroduced in the “State of the Art”, and is referred to as GMUSIC. Asopposed to the traditional MUSIC algorithm, GMUSIC is designed toprovide consistent estimates of the DoAs when both the sample size andthe number of sensors or antennas have the same order of magnitude. Thisensures a significant improvement of the performance in the low powerand/or low sample size regime.

The MUSIC algorithm provides consistent results of the estimateddirections-of-arrival or angles-of-arrival when N is high but M remainsfixed, that is to say, when N is high relative to M. On the contrary,our algorithm provides consistent and reliable results when both N and Mare high, that is to say, when M and N are comparable to each other interms of order of magnitude. In practice, this means that our algorithmbehaves in a very reliable way even if the number of samples N is notmuch higher than the number of sensors or antennae M. This does nothappen in the original MUSIC algorithm. Therefore, our algorithm behavesmuch better than the MUSIC algorithm in the “low-sample size regime”.

The same reasoning applies to the low-power regime. It is well knownthat, in every algorithm of directions-of-arrival detection or spectralestimation, the higher the power of the signals to be detected is, theless number of samples is needed in order to have a reliable estimationof the DoA of the signals. Since our algorithm allows good detectionwith a low number of samples with respect to the MUSIC algorithm, we areable to detect signals whose power is low with respect to the signalpower which the MUSIC algorithm would require in order to detect thosesignals, for the same number of samples.

The method comprises the following steps:

-   -   1. From the collection of N column vectors or snapshots, {y(1),        . . . , y(N)}, calculate a sample spatial correlation matrix or        sample covariance matrix {circumflex over (R)}.        -   In a preferred embodiment, this sample spatial correlation            matrix takes the following form:

$\begin{matrix}{\hat{R} = {\frac{1}{N}{\sum\limits_{n = 1}^{N}{{y(n)}{y^{H}(n)}}}}} & (2)\end{matrix}$

-   -   -    where (•)^(H) denotes transpose if we are dealing with            real-valued signals, and transpose conjugate if we are            dealing with complex-valued signals. However, this sample            spatial correlation matrix is not limited to this form (2),            but the sample spatial correlation matrix {circumflex over            (R)} can take any other form from the ones known in the art.

    -   2. Obtain the eigenvalues and eigenvectors of the matrix        {circumflex over (R)}. They are denoted as {ê_(k), k=1 . . . M}        and {circumflex over (λ)}₁≦{circumflex over (λ)}₂≦ . . .        ≦{circumflex over (λ)}_(M), respectively.

    -   3. Obtain all the distinct solutions to the following equation        in μ:

$\begin{matrix}{{\sum\limits_{k = 1}^{M}\frac{{\hat{\lambda}}_{k}}{{\hat{\lambda}}_{k} - \mu}} = N} & (3)\end{matrix}$

-   -    these solutions are ordered as μ₁≦ . . . ≦μ_(M).        -   If M≦N, these M solutions are, in principle, all different.            If some of values of {circumflex over (λ)}_(k) are not            different, that is to say, they are repeated in the list            {{circumflex over (λ)}₁, . . . , {circumflex over (λ)}_(M)},            then the value μ_(k) is also repeated in the list {μ₁, . . .            , μ_(M)} as many times as {circumflex over (λ)}_(k) is also            repeated in the list {{circumflex over (λ)}₁, . . . ,            {circumflex over (λ)}_(M)}.        -   If M>N, we set μ₁= . . . =μ_(M−N)=0, so that μ_(M−N+1)≦ . .            . ≦μ_(M) contain the positive solutions to said equation in            μ. These solutions are, in principle, all different. Again,            if some of values of {circumflex over (μ)}_(k) are not            different, that is to say, they are repeated in the list            {{circumflex over (μ)}_(M−N+1), . . . , {circumflex over            (μ)}_(M)}, then the value μ_(k) is also repeated in the list            {μ_(M−N+1), . . . , μ_(M)} as many times as {circumflex over            (λ)}_(k) is also repeated in the list {{circumflex over            (λ)}_(M−N+1), . . . , {circumflex over (λ)}_(M)}.    -   4. Construct the following weights

$\begin{matrix}{\varphi_{m} = \left\{ \begin{matrix}1 & {m \leq \left\lbrack {M - N} \right\rbrack^{+}} \\{1 - \Phi_{m} + \Psi_{m}} & {{\left\lbrack {M - N} \right\rbrack^{+} + 1} \leq m \leq {M - K}} \\{- \Phi_{m}} & {m > {M - K}}\end{matrix} \right.} & (4) \\{\Phi_{m} = {\sum\limits_{k = {{\lbrack{M - N}\rbrack}^{+} + 1}}^{M - K}\frac{{\hat{\lambda}}_{m}\left( {{\hat{\lambda}}_{k} - \mu_{k}} \right)}{\left( {{\hat{\lambda}}_{m} - {\hat{\lambda}}_{k}} \right)\left( {{\hat{\lambda}}_{m} - \mu_{k}} \right)}}} & (5) \\{{\Psi_{m} = {\sum\limits_{\underset{k \neq m}{k = {{\lbrack{M - N}\rbrack}^{+} + 1}}}^{M}\frac{{\hat{\lambda}}_{k}\left( {{\hat{\lambda}}_{m} - \mu_{m}} \right)}{\left( {{\hat{\lambda}}_{k} - {\hat{\lambda}}_{m}} \right)\left( {{\hat{\lambda}}_{k} - \mu_{m}} \right)}}}{where}{{{M - N}}^{+} = \left\{ {\begin{matrix}{M - N} & {{{if}\mspace{14mu} M} > N} \\0 & {{otherwise}.}\end{matrix}.} \right.}} & (6)\end{matrix}$

-   -   5. The estimated DoAs are chosen as either by selecting the K        deepest local minima of the cost function ƒ(η(θ)), where ƒ(•) is        an increasing real-valued function, or by selecting the K        highest local maxima of the following cost function g(η(θ)),        where g(•) is a decreasing real-valued function, and wherein

$\begin{matrix}{{\eta(\theta)} = {{a^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){{a(\theta)}.}}} & (7)\end{matrix}$

Observe the similarity with the traditional MUSIC method, which selectsthe estimated DoAs as the K deepest local minima of the function

$\begin{matrix}{{\eta_{MUSIC}(\theta)} = {{a^{H}(\theta)}\left( {\sum\limits_{k = 1}^{M - K}{{\hat{e}}_{k}{\hat{e}}_{k}^{H}}} \right){{a(\theta)}.}}} & (8)\end{matrix}$

The present invention exploits the information from all theeigenvectors, not only of those associated with the lowest eigenvalues(as it is the case in the traditional MUSIC method). Next, it will beshown that the present invention overcomes the performance breakdown ofthe MUSIC method in low power or low sample size scenarios.

FIG. 3 is a graph showing a particular instance of the functions in (7)and (8) above as a function of the angle of arrival, in a situationwhere two point sources were located at an azimuth of 35° and 40° andtheir signals were received in Gaussian white noise. The sourcesilluminated a uniform linear array of sensors or antennas with M=10elements disposed linearly with an inter-element separation equal tohalf the wavelength of the impinging waveforms. In this example, a totalof N=70 snapshots were available at the receiver, and the signals werereceived with a power 10 dB higher than the background noise power. Inthis particular realisation of the experiment, the traditional MUSICmethod is unable to distinguish the two directional sources, and thedeepest local minima were located at 39.51° and −32.82°. With theproposed GMUSIC invention, the system is able to recognize two angularsources at 34.92° and 39.58°, very close to the actual values of thesource locations.

FIG. 4 compares the mean squared error obtained by averaging 500realisations of the proposed and traditional MUSIC methods when N=7 andN=70 snapshots were available at the array. The scenario was the same asin FIG. 3. Results are shown for different value of the signal-to-noiseratio (SNR), which is the quotient between the power of the receivedsignals (assumed equal for the two sources) and the power of thebackground omni directional noise. The unconditional Cramér-Rao Bound(CRB), that is to say, the lowest variance that can be achieved with anunbiased estimator of the direction of arrival, is also shown forinformative purposes. Observe that GMUSIC method presents an SNRthreshold much lower than the traditional MUSIC. Furthermore, theproposed approach shows also a better performance in the high SNRregion.

FIG. 5 represents the probability of detection of the traditional MUSICand the GMUSIC systems, taken from an average of 10⁴ realisations. Theprobability of detection defined as

${{Prob}\left\lbrack {{\eta\left( \frac{\theta_{1} + \theta_{2}}{2} \right)} > \frac{{\eta\left( \theta_{1} \right)} + {\eta\left( \theta_{2} \right)}}{2}} \right\rbrack},$where θ₁ and θ₂ are the DoAs of two point sources in the simulatedscenario. In particular, we fixed θ₁=0° whereas θ₂ varied from 0° to10°. Results are shown for two values of the sample size, N=7 and N=70,and two values of the SNR, namely SNR=10 dB and SNR=40 dB. It isexperimentally demonstrated that the GMUSIC method also overcomes thetraditional MUSIC in terms of probability of detection of two closelyspaced sources.

As already mentioned in this description, our method can be furthercomplemented with some enhancements of the traditional MUSIC technique,such as polynomial rooting, beamspace projection, sequential estimationor spatial smoothing, which are shown in the following preferredembodiments.

Thus, when the method is applied to a uniform linear array or receivingelements, the DoAs can also be estimated using a trivial modification ofour method as previously defined, as follows. This procedure is usuallyreferred to as polynomial rooting.

Define the M×M matrix

$\Gamma = {\sum\limits_{k = 1}^{M}{\varphi_{k}{\hat{e}}_{k}{\hat{e}}_{k}^{H}}}$and let Γ_(m,n) denote the entry corresponding to the mth row and thenth column.

Consider the complex-valued polynomial

${D(z)} = {\sum\limits_{l = {- {({M - 1})}}}^{M - 1}{\alpha_{l}z^{M - l - 1}}}$wherein $\alpha_{l} = \left\{ \begin{matrix}{\sum\limits_{m = 1}^{M - l}\Gamma_{m,{m + l}}} & {l = {0\mspace{14mu}\ldots\mspace{14mu}\left( {M - 1} \right)}} \\\alpha_{- l}^{*} & {l = {{{- \left( {M - 1} \right)}\mspace{14mu}\ldots}\; - 1}}\end{matrix} \right.$and where (•)* denotes complex conjugate.

The DoAs are estimated as the values

$\theta_{1} = {{arc}\;{\sin\left( {\frac{\lambda}{2\;\pi\; d}{\arg\left( z_{1} \right)}} \right)}}$⋮$\theta_{K} = {{arc}\;{\sin\left( {\frac{\lambda}{2\;\pi\; d}{\arg\left( z_{k} \right)}} \right)}}$where d denotes the distance between the elements of the array, λ is thewavelength of the impinging waveforms, arg(•) denotes the argument of acomplex number, and where z₁ . . . z_(K) are chosen out of the 2(M−1)roots of the complex polynomial D(z), denoted by {z(k), k=1 . . .2(M−1)}, as those that either minimize the quantityƒ(ƒD(exp(jarg(z_(k))))|), wherein ƒ(•) is an increasing real-valuedfunction, or maximize the quantity g(|D(exp(jarg(z_(k))))|), whereing(•) a decreasing real-valued function, and wherein j is the imaginaryunit, arg(•) denotes the argument of a complex number and ∥ denotes themodule of a complex number.

Another procedure that can be incorporated into the method as previouslyproposed herein, when applied to uniform linear arrays of receivingelements, is spatial smoothing. This method consists in devoting spatialdegrees of freedom to improving the estimation of the spatial covariancematrix. In this technique, the array of M elements is divided intoM_(ss) smaller sub-arrays, each containing L sensors or antennas. Ifthere is no overlapping between the M_(ss) smaller sub-arrays, thenM_(ss)×L=M. If, on the contrary, there is overlapping between the M_(ss)smaller sub-arrays, then M_(ss)×L>M. Each of the N column vectors y(1),. . . , y(N) is divided into M_(ss) smaller column vectors y_(Mss)(1), .. . , y_(Mss)(N), each of said smaller column vectors having L rows;

The covariance matrix is estimated as follows. First, the covariance orsample spatial correlation matrix at each of the M_(ss) subarrays isestimated as the corresponding sample correlation matrix in (2) or usingany other estimation of the spatial correlation matrix as known in theart. That is to say, M_(ss) sample spatial correlation matrices at eachof the M_(ss) smaller subarrays are calculated. Second, an L×Lcovariance matrix estimate {circumflex over (R)} is constructed byaveraging the M_(ss) covariance matrices obtained in the previous step.Let α_(ss)(θ) denote a column vector formed by selecting the first Lrows of α(θ). The method consists of replacing α(θ) with α_(ss)(θ) andthe sample covariance matrix {circumflex over (R)} by the estimationdescribed above. Then, either K deepest local minima of the followingcost function ƒ(η(θ)), wherein f(•) is an increasing real-valuedfunction, or K highest local maxima of the following cost functiong(η(θ)), where g(•) is a decreasing real-valued function, arecalculated, wherein

${\eta(\theta)} = {{a_{ss}^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){{a_{ss}(\theta)}.}}$

Another possible modification of our algorithm as previously described,that may lead to further improved results, follows from the sequentialdetection approach. Instead of selecting the K deepest local minima (orK highest local maxima) of the cost function in (7), namely

${{\eta(\theta)} = {{a^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){a(\theta)}}},$the approach estimates the DoAs following a process of K iterations,wherein at the kth iteration, k=1 . . . K, the kth DoA is selected aseither the corresponding minimum of the cost function ƒ(ν_(k)(θ)),wherein where f(•) is an increasing real-valued function, or thecorresponding maximum of the cost function g(ν_(k)(θ)), wherein g(•) isa decreasing real-valued function, and wherein

${\vartheta_{k}(\theta)} = {{b_{k}^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){b_{k}(\theta)}}$where${{b_{k}(\theta)} = \frac{\left\lbrack {I - {{A_{k}\left( {A_{k}^{H}A_{k}} \right)}^{- 1}A_{k}^{H}}} \right\rbrack{a(\theta)}}{{\left\lbrack {I - {{A_{k}\left( {A_{k}^{H}A_{k}} \right)}^{- 1}A_{k}^{H}}} \right\rbrack{a(\theta)}}}},{A_{k} = \left\lbrack {{a\left( \theta_{1} \right)},\ldots\mspace{14mu},{a\left( \theta_{k - 1} \right)}} \right\rbrack}$and where θ₁, . . . , θ_(k−1) are the DoAs detected in previousiterations and ∥•∥ denotes the Euclidean norm of a vector.

The method proposed herein is also applicable to the case where a lineartransformation is applied to the observations or the sample covariancematrix in order to improve the properties of the DoA estimationprocedure.

For example, it is common practice in this type of applications tomultiply the receiving snapshots by a M_(b)×M matrix B, where M<M_(b).This procedure is usually referred to as beamspace projection. Thematrix B is designed taking into account the array geometry in order toenhance the resolution of the algorithm across certain regions of DoAsof interest. In other words, the values of the matrix B depend on thearray geometry of the receiving elements. Our method can be triviallymodified to operate in this situation. Instead of using the samplecovariance matrix {circumflex over (R)} in (2), one must use thefollowing M_(b)×M_(b) transformed version thereof

${\hat{R}}_{B} = {{B^{H}\left( {\frac{1}{N}{\sum\limits_{n = 1}^{N}{{y(n)}{y^{H}(n)}}}} \right)}B}$

From this point, the DoAs of the signals of interest are estimated byeither selecting K deepest local minima of the cost functionƒ(η_(b)(θ)), wherein f(•) is an increasing real-valued function, orselecting K highest local maxima of the cost function g(η_(B)(θ)),wherein g(•) is a decreasing real-valued function, and wherein

${\eta_{B}(\theta)} = \frac{{a^{H}(\theta)}{B\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}^{B}{{\hat{e}}_{m}^{B}\left( {\hat{e}}_{m}^{B} \right)}^{H}}} \right)}B^{H}{a(\theta)}}{{a^{H}(\theta)}{BB}^{H}{a(\theta)}}$where {ê_(m) ^(B), m=1 . . . M} denote the eigenvectors of {circumflexover (R)}_(B) and where the weights {φ_(m) ^(B), m=1 . . . M} areconstructed as φ_(m) replacing {circumflex over (R)} with {circumflexover (R)}_(B).

In this text, the term “comprises” and its derivations (such as“comprising”, etc.) should not be understood in an excluding sense, thatis, these terms should not be interpreted as excluding the possibilitythat what is described and defined may include further elements, steps,etc.

The invention is obviously not limited to the specific embodimentsdescribed herein, but also encompasses any variations that may beconsidered by any person skilled in the art (for example, as regards thechoice of components, configuration, etc.), within the general scope ofthe invention as defined in the appended claims.

1. A method for determining the direction of arrival of at least oneradio waveform (w_(a)), comprising the steps of: receiving at least onewaveform (w_(a)) by M receiving elements (1, 2, 3, 4 . . . ), M being anatural number and M>1, and obtaining one detected signal (d_(i), d₂,d₃, d₄ . . . ) at the output of each of the M receiving elements (1, 2,3, 4 . . . ); processing said M detected signals (d_(i), d₂, d₃, d₄ . .. ) at M pre-processing units (30, 31, 32, 33 . . . ) for obtaining Mprocessed signals (f₁, f₂, f₃, f₄ . . . ); sampling at N time instantssaid M processed signals (f₁, f₂, f₃, f₄ . . . ) at a signal processingunit (40), N being a natural number, for obtaining corresponding samplevalues; from the sample values obtained in the previous step, creating Ncolumn vectors y(1), . . . , y(N), wherein y(n) represents, for asampling instant n, a column vector with M rows, each of the M rowscontaining the sample value associated to the corresponding receivingelement (1, 2, 3, 4 . . . ) at that sampling instant n; from thecollection of N column vectors y(1), . . . , y(N), calculating a samplespatial correlation matrix {circumflex over (R)}; obtaining all theeigenvectors ê_(k), where k=1 . . . M, and all the eigenvalues{circumflex over (λ)}₁≦{circumflex over (λ)}₂≦ . . . ≦{circumflex over(λ)}_(M) from said sample spatial correlation matrix {circumflex over(R)}; the method being characterised in that it further comprises thesteps of: from all the eigenvectors ê_(k) and all the eigenvalues{circumflex over (λ)}₁≦{circumflex over (λ)}₂≦ . . . ≦{circumflex over(λ)}_(M) obtained in the previous step, obtaining a set of signaldependent parameters μ_(k), where k=1 . . . M; from the eigenvalues{circumflex over (λ)}₁≦{circumflex over (λ)}₂≦ . . . ≦{circumflex over(λ)}_(M) and the signal dependent parameters μ_(k) previously obtained,constructing a set of weights; finding K solutions from a certain costfunction such that said K solutions fulfil a certain requirement.
 2. Themethod according to claim 1, wherein N>M, N and M being natural numbers.3. The method according to claim 1, wherein N=M, N and M being naturalnumbers.
 4. The method according to claim 1, wherein N<M, N and M beingnatural numbers.
 5. The method according to claim 1, wherein the step ofreceiving at least one waveform (w_(a)) is done by an array of M sensoror antenna elements (1, 2, 3, 4 . . . ).
 6. The method according toclaim 1, wherein there are K sources, each of them emitting a waveform.7. The method according to claim 1, further comprising the step ofacquiring and storing, at the signal processing unit (40), the values ofa column vector a(θ), said column vector a(θ) having M rowscorresponding to respective M receiving elements, each of the M rowscomprising a function which depends on the direction of arrival θ thatis to be explored.
 8. The method according to claim 7, furthercomprising the step of calibrating the array of M receiving elements asfollows: exposing a collection of reference signals to the array of Mreceiving elements, said reference signals simulating waveforms comingfrom the direction θ; measuring, at the signal processing unit (40), theresponse of the received reference signals obtained at each of the Mreceiving elements, thereby obtaining a(θ).
 9. The method according toclaim 1, further comprising the step of calibrating the array of Mreceiving elements (1, 2, 3, 4 . . . ) as follows: exposing a collectionof reference signals to the array of M receiving elements (1, 2, 3, 4 .. . ); measuring, at the signal processing unit (40), the response ofthe received reference signals obtained at the output of the Mpre-processing units (30, 31, 32, 33, . . . ); modifying the response ofthe M pre-processing units (30, 31, 32, 33, . . . ) in such a way thatthe signals measured at the output of the M pre-processing units (30,31, 32, 33, . . . ) are as close as possible to the reference signalsinserted in the first step.
 10. The method according to claim 1, whereinthe step of obtaining a set of signal dependent parameters μ_(k) is doneby obtaining all the distinct solutions μ₁≦ . . . ≦μ_(M) to thefollowing equation in μ:${\sum\limits_{k = 1}^{M}\frac{{\hat{\lambda}}_{k}}{{\hat{\lambda}}_{k} - \mu}} = {N.}$11. The method according to claim 10, wherein if N<M, μ₁= . . .=μ_(M−N)=0 and wherein μ_(M−N+1)≦ . . . ≦μ_(M) are the positivesolutions to said equation in μ.
 12. The method according to claim 1,wherein the step of constructing a set of weights φ_(m) is done byapplying the following equation: $\varphi_{m} = \left\{ \begin{matrix}1 & {m \leq \left\lbrack {M - N} \right\rbrack^{+}} \\{1 - \Phi_{m} + \Psi_{m}} & {{\left\lbrack {M - N} \right\rbrack^{+} + 1} \leq m \leq {M - K}} \\{- \Phi_{m}} & {m > {M - K}}\end{matrix} \right.$ wherein K is the number of sources and thereforethe number of waveforms (w_(a), w_(b), w_(c) . . . ) received at the Mreceiving units (1, 2, 3, 4 . . . ) and wherein:$\Phi_{m} = {\sum\limits_{k = {{\lbrack{M - N}\rbrack}^{+} + 1}}^{M - K}\frac{{\hat{\lambda}}_{m}\left( {{\hat{\lambda}}_{k} - \mu_{k}} \right)}{\left( {{\hat{\lambda}}_{m} - {\hat{\lambda}}_{k}} \right)\left( {{\hat{\lambda}}_{m} - \mu_{k}} \right)}}$$\Psi_{m} = {\sum\limits_{\underset{k \neq m}{k = {{\lbrack{M - N}\rbrack}^{+} + 1}}}^{M}\frac{{\hat{\lambda}}_{k}\left( {{\hat{\lambda}}_{m} - \mu_{m}} \right)}{\left( {{\hat{\lambda}}_{k} - {\hat{\lambda}}_{m}} \right)\left( {{\hat{\lambda}}_{k} - \mu_{m}} \right)}}$${{where}\left\lbrack {M - N} \right\rbrack}^{+} = \left\{ \begin{matrix}{M - N} & {{{if}\mspace{14mu} M} > N} \\0 & {{otherwise}.}\end{matrix} \right.$
 13. The method according to claim 1, wherein thesample spatial correlation matrix {circumflex over (R)} takes the formof:$\hat{R} = {\frac{1}{N}{\sum\limits_{n = 1}^{N}{{y(n)}{{y^{H}(n)}.}}}}$14. The method according to claim 1, wherein the sample spatialcorrelation matrix {circumflex over (R)} takes any other form from theones known in the art.
 15. The method according to claim 12, wherein thestep of finding the K solutions from the cost function is done either byselecting K deepest local minima of the following cost function:f(η(θ)) where f(•) is an increasing real-valued function, or byselecting the K highest local maxima of the following cost function:g(η(θ)) where g(•) is a decreasing real-valued function, and wherein${{\eta(\theta)} = {{a^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){a(\theta)}}};{and}$further comprising the step of acquiring and storing, at the signalprocessing unit (40), the values of a column vector a(θ), said columnvector a(θ) having M rows corresponding to respective M receivingelements, each of the M rows comprising a function which depends on thedirection of arrival θ that is to be explored.
 16. The method accordingto claim 12 wherein, when the M receiving elements (1, 2, 3, 4 . . . )form a uniform linear array, the method further comprises the steps of:constructing the matrix${\Gamma = {\sum\limits_{k = 1}^{M}{\varphi_{k}{\hat{e}}_{k}{\hat{e}}_{k}^{H}}}};$constructing a polynomial${D(z)} = {\sum\limits_{l = {- {({M - 1})}}}^{M - 1}{\alpha_{l}z^{M - l - 1}}}$wherein $\alpha_{l} = \left\{ \begin{matrix}{\sum\limits_{m = 1}^{M - l}\Gamma_{m,{m + l}}} & {l = {0\mspace{14mu}\ldots\mspace{14mu}\left( {M - 1} \right)}} \\\alpha_{- l}^{*} & {{l = {{{- \left( {M - 1} \right)}\mspace{14mu}\ldots} - 1}},}\end{matrix} \right.$ and (•)* denotes complex conjugate and Γ_(m,n)denotes the entry corresponding to the mth row and the nth column of thematrix Γ; choosing z₁, . . . , z_(K) from the roots of the complexpolynomial D(z), as those that either minimize the quantityf(|D(exp(jarg(z_(k))))|), wherein f(•) is an increasing real-valuedfunction, or maximize the quantity g(|(exp(jarg(z_(k))))|) wherein g(•)a decreasing real-valued function, and wherein j is the imaginary unit,arg(•) denotes the argument of a complex number and ∥ denotes the moduleof a complex number; obtaining, from those z₁, . . . , z_(K), the Kestimated directions-of-arrival:$\theta_{1} = {\arcsin\left( {\frac{\lambda}{2\pi\; d}\mspace{14mu}{\arg\left( z_{1} \right)}} \right)}$⋮$\theta_{K} = {\arcsin\left( {\frac{\lambda}{2\pi\; d}\mspace{14mu}{\arg\left( z_{K} \right)}} \right)}$wherein d denotes the distance between the elements of the array and λis the wavelength of the impinging waveforms.
 17. The method accordingto claim 12, wherein, when the M receiving elements (1, 2, 3, 4 . . . )form a uniform linear array, the method further comprises the steps of:dividing the array of M elements into M_(ss) smaller subarrays, eachcontaining L receiving elements, M_(ss) and L being natural numbers;dividing each of the N column vectors y(1), . . . , y(N) into M_(ss)smaller column vectors Y_(Mss)(1), . . . , Y_(Mss)(N), each of saidsmaller column vectors having L rows; calculating M_(ss) sample spatialcorrelation matrices at each of the M_(ss) smaller subarrays; averagingsaid M_(ss) sample spatial correlation matrices and constructing anaveraged L×L spatial correlation matrix {circumflex over (R)} therefrom;constructing a column vector a_(ss)(θ) by selecting the first L rows ofa(θ); selecting either K deepest local minima of the following costfunction:f(η(θ)) wherein f(•) is an increasing real-valued function, or K highestlocal maxima of the following cost functiong(η(θ) where g(•) is a decreasing real-valued function, and${{\eta(\theta)} = {{a^{H}(\theta)}\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{\hat{e}}_{m}{\hat{e}}_{m}^{H}}} \right){a(\theta)}}};{and}$further comprising the step of acquiring and storing, at the signalprocessing unit (40), the values of a column vector a(θ), said columnvector a(θ) having M rows corresponding to respective M receivingelements, each of the M rows comprising a function which depends on thedirection of arrival θ that is to be explored.
 18. The method accordingto claim 12, wherein the step of finding K solutions from a certain costfunction is done by selecting K local extrema of said cost function bymeans of a process of K iterations, wherein at the kth iteration, k=1 .. . K, the kth local extremum is selected as either the correspondingminimum of the cost functionf(θ_(k)(θ)) wherein where f(•) is an increasing real-valued function, orthe corresponding maximum of the cost functiong(θ_(k)(θ)) wherein g(•) is a decreasing real-valued function, andwherein${{b_{k}(\theta)} = \frac{\left\lbrack {I - {{A_{k}\left( {A_{k}^{H}A_{k}} \right)}^{- 1}A_{k}^{H}}} \right\rbrack{a(\theta)}}{{\left\lbrack {I - {{A_{k}\left( {A_{k}^{H}A_{k}} \right)}^{- 1}A_{k}^{H}}} \right\rbrack{a(\theta)}}}},{A_{k}\left\lbrack {{a\left( \theta_{1} \right)},\ldots\mspace{11mu},{a\left( \theta_{k - 1} \right)}} \right\rbrack}$and wherein θ₁, . . . , θ_(k−1) are the local extrema selected inprevious iterations and ∥•∥ denotes the Euclidean norm of a vector; andfurther comprising the step of acquiring and storing, at the signalprocessing unit (40), the values of a column vector a(θ), said columnvector a(θ) having M rows corresponding to respective M receivingelements, each of the M rows comprising a function which depends on thedirection of arrival θ that is to be explored.
 19. The method accordingto claim 12, further comprising the step of enhancing the resolution ofthe algorithm across selected regions of directions-of-arrival, bymultiplying the sample values obtained at the signal processing unit(40) by a M_(b)×M matrix B, where M<M_(b), wherein the values of thematrix B depend on the array geometry of the receiving elements.
 20. Themethod according to claim 19, wherein the sample spatial correlationmatrix {circumflex over (R)} takes the form of:${\hat{R}}_{B} = {{B^{H}\left( {\frac{1}{N}{\sum\limits_{n = 1}^{N}{{y(n)}{y^{H}(n)}}}} \right)}{B.}}$21. The method according to claim 20, wherein the directions-of-arrivalof the signals of interest are estimated by either selecting K deepestlocal minima of the cost functionf(η_(B)(θ)) wherein f(•) is an increasing real-valued function, orselecting K highest local maxima of the cost functiong(η_(B)(θ)) wherein g(•) is a decreasing real-valued function, andwherein${\eta_{B}(\theta)} = {\frac{{a^{H}(\theta)}{B\left( {\sum\limits_{m = 1}^{M}{\varphi_{m}{{\hat{e}}_{m}\left( {\hat{e}}_{m} \right)}^{H}}} \right)}B^{H}{a(\theta)}}{{a^{H}(\theta)}{BB}^{H}{a(\theta)}}.}$22. The method according to claim 1, wherein said at least one receivingwaveform (w_(a)) is a radio-frequency signal or a microwave signal. 23.A system comprising means adapted for carrying out the steps of themethod according to claim
 1. 24. The system according to claim 23,wherein the receiving units (1, 2, 3, 4 . . . ) are an array of M sensoror antenna elements.
 25. The system according to claim 23, wherein the Mpre-processing units (30, 31, 32, 33 . . . ) are radio-frequency oracoustic front-ends (30, 31, 32, 33 . . . ).
 26. The system according toclaim 23, wherein the at least one received waveform (w_(a)) isgenerated by at least one source.
 27. The system according to claim 26,wherein said at least one source is a transmitter.
 28. The systemaccording to claim 26, wherein said at least one source is a reflectorof the acoustic or radio waveform generated by other sources.
 29. Thesystem according to claim 23, wherein said at least one receivingwaveform (w_(a)) is a radio-frequency signal or a microwave signal. 30.A computer program comprising computer program code means adapted toperform the steps of the method of claim 1 when said program is run on acomputer, a digital signal processor, a field-programmable gate array,an application-specific integrated circuit, a micro-processor, amicro-controller, or any other form of programmable hardware.